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ABSTRACT 

Real world phenomena — exhibit nonlinear relationships, complex 
geometry, and intricate processes. Analytic or exact solution methods only address a 
minor class of such phenomena. Consequently, numerical approximation methods, such 
as pootiading methods, can be used. 

The goal is, by making use of a variety of root-finding methods Keven 
Rhapson, Chebyshev, Halley and iaeuerie): to gain a qualitative appreciation on how 
various root-finding methods address many prevailing real-world concerns, to include, 
how are suitable approximation methods determined; when do root finding methods 
converge; and how long for convergence? 

Answers to the questions were gained through examining the basins of attraction 
of the root-finding methods. Different methods generate different basins of attraction. In 


the end, each method appears to have its own advantages and disadvantages. 
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I. INTRODUCTION 


A. BACKGROUND 

Root finding methods have been of interest for a long time. Why? Often people 
ask qualitative questions about real-world phenomena, and they want these questions 
answered. To come to an answer, one must accurately model the real world phenomena 
in a mathematical model, and then solve the model. In many applications, the solution 
involves finding a root. 

Constructing models is rarely a simple process. Models come in many shapes and 
sizes. Some of these represent a dynamical process — a recipe for how real-world 
phenomena interact and change over time. How these interactions and changes occur 
governs the choice of model. For example, the continuous model leading to a differential 
equation is reasonable for certain phenomena, while difference equations in the form of a 
recurrence relation address phenomena occurring in discrete steps. Solutions, however, 
are not guaranteed in every instance. 

When analytical or exact methods are applicable, sometimes formulas for 
solutions exist. However, these methods are restrictive, often providing insight into the 
behavior of only a minor class of real world phenomena. Included in this category are 
models that can be approximated by linear relationships, simple geometry, and low 
dimensionality. For a great deal of real world phenomena, that is not the norm. Real 
ai phenomena commonly exhibit nonlinear relationships, complex geometry, and 
intricate processes. Consequently, exact methods can be of limited practical value 


(Chapra, 1988). 





B. MOTIVATION 

Where analytical or exact methods fail, numerical approximation methods often 
succeed, approximately. One such approximation method employs difference equations. 
When applied to a large though finite number of steps, difference equations are closely 
related to the continuous behavior of a differential equation (Figure 1.1) In fact, a 


continuous model, y(t), can be seen as a limit of the discrete model, y,(¢,) (Figure 1.2). 


J / | 
Xp Xy X> X3 Xo X X> X3 Xq Xs X6 Xo aoe een ere eee aeres X, 


Figure 1.1. Approximating Continuous Behavior 


2 = f(t,y) 


_ lim Vnel mee 
Aat>0 At 


mw Yast Yn 
At 


Dra eames 
ar AO mal 
F(t, +¥,) ) 


Y nel Vin + Atl f(t,.y, JI 


Figure 1.2. Approximating a Differential Equation Using a Difference Equation 
Although the model approaches are different, solution methods for each share 


common ground. In the contmuous model, solution curves may be obtained from the 











roots of a linear, constant-coefficient differential equation’s characteristic polynomial. In 
the discrete walik solutions come from the roots of the recurrence relation’s 
characteristic polynomial. In either case, roots can be real, imaginary, or complex. 
Consequently, solutions can vary greatly in their dynamical behaviors. 

Numerical (Root finding) methods, however, serve as the computational tools that 
unveil the mysteries of such dynamical behavior. Different methods, however, may 
produce different results from the same initial guess. So things can get really interesting! 
C. GOALS 

This thesis seeks to gain a qualitative appreciation on how various root-finding 
methods address many prevailing real-world concerns, to include, how are suitable 
approximation methods determined; when do root finding methods converge; and how 
long for convergence? Particular emphasis is given to finding which initial guesses lead 
to which roots. 

D. METHODOLGY 

From a mesh of points within the complex plane, Newton-Raphson, Chebyshev, 
Halley, and Laguerre root-finding methods numerically compute the successive 
approximations of some n” order complex polynomial’s roots. In order to better grasp 
the effects, the results are mapped — thus, creating a geometry of basins of attractions 
which are the set of starting points whose trajectories are asymptotic to a bounded region 
— (Devaney, 1989). The geometrical differences lend to the qualitative difference amongst 


the root-finding methods. 
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Il. BEHAVIORS OF DYNAMICAL SYSTEMS — 


While the mathematics describing dynamic behavior may be fairly 
straightforward, interpreting such behavior can be difficult. In order to truly grasp it, one 
must familiarize oneself with the role of numerical methods and the utility of mapping 
their geometry. | 
A. NUMERICAL METHODS ROLE 

Numerical methods approximate solutions to mathematically expressed models. 
When rane solutions are obtained from the zeros of some functions, root-finding 
methods serve as the tool of choice. These methods are usually iterative — beginning with 
an initial starting value and computing successive approximate solutions using a well- 


defined recurrence relation (Figure 2.1). Each successive step yields a numerical solution 


Sample Recurrence Relation Idea of Successive Approximations 
Xn = F(%) Xx,=F CF Co Fl%) -) 


n times 


Figure 2.1. Recurrence Relation & Iteration 


to the recurrence relation — in essence, generating a sequence of even better 
approximations. Hence, the solution process itself is a discrete dynamical system that 


generates a sequence of numbers. As Table 2.2 illustrates, each term of the sequence not 


Sequence Numerical Solutions Per Iterative Step Long Term 
Behavior 


J | ee ee (a en ey ee Q (Convergence) 
H Sy ee ey A Sk ak i co (Divergence) 
Ip Ps-VarTarVarVatrarn >? 
Figure 2.2. Arbitrary Numerical Solutions 
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only signifies a numerical clades for the n” iterative step, but also suggests the 
ultimate behavior of that solution. Determining such behavior is not always done by 
simple inspection. Some sequences are obvious; others are not. Consequently, numbers 
alone are often not enough. 

B. UTILITY OF MAPPING — SINGLE FIXED POINT 

Another, often preferred, method used to determine dynamical behaviors is to 
visualize them. Visualization entails mapping out the geometry of the numerical 
solutions. Why is this geometry important? Simply put, it graphically depicts the 
dynamical behavior of root-finding methods. 

As Figure 2.3 suggests, the mapping of a sequence of numerical solutions depicts 
the behavioral path or trajectory of a single starting point. Starting points that do not 
change after iteration are called fixed, and qualitative behaviors of other starting points 
can be interpreted in relation to the fixed points. 


Sequence I Sequence II _ Sequence III 











Figure 2.3. Mapping of Numerical Solutions 


Cobweb diagrams help point out the qualitative behaviors near fixed points using 
the principle of feedback (Figure 2.4) (After Peitgen, 1992). The principle of feedback is 
simple — an input, x,, is given, processed through some function, f, and then the 


output, y,, becomes the next input, x,+;, repeatedly. When allowing the ouput to equal 
6 





the next input, an identity exists so that (x,,, = y,) =(x = y). Cobweb diagrams exploit 


+] 


the relationship, map the iterations, and reveal the behaviors of fixed points. 





Figure 2.4. Cobweb Diagram from the Principle of Feedback 


Behaviors about fixed points are converging, diverging or chaotic, and all can be 
mapped. Convergent mappings point out attracting fixed points; divergent mappings 
denote repulsive ones; and chaotic mappings never settle. While all three behaviors are 
essential in describing dynamical behavior, a simple example excluding chaotic 


mappings will suffice to illustrate the point. Consider a simple linear recurrence 


<1, the mapping contracts, 








relation; x,,,= f(x,)=mx,+b. When 


I) 





converging to a fixed point. When |/ (x,)| > 1, the mapping expands diverging off to 


positive or negative infinity. Figure 2.5 clarifies the point. 

With relatively little effort, the geometrical approach can handle nonlinear 
behavior as well. For smooth nonlinear recurrence relations, the same sort of contracting 
ee argument holds near the fixed point. As Figure 2.6 indicates, the trick is 
to locally reduce the nonlinear model to linear parts, apply the graphical analysis, and 


then couple the pieces together. 
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Figure 2.6. Concept of Linearizing a Nonlinear Mapping Function 








Although dynamic behaviors about a single fixed point are fairly predictable, 
startling behavioral effects can and often do occur when multiple fixed points exist. 
C. UTILITY OF MAPPING — MULTIPLE FIXED POINT 

When fixed points coexist, the geometry of numerical solutions can change 
considerably. The effect of each fixed point is no longer simple; rather their effects 
interact. Consequently, determining such points and their effect is a necessity. 

Multiple fixed points are often found in the realm of nonlinear phenomena. While - 
the effect of a single fixed point has been discussed, what happens when there are two, 
three, or n of them? How many can there be when the iterator is a root approximation 
method? As Figure 2.7 points out, the fundamental theorem of algebra tells us that an n” 
degree polynomial is factorable into n linear factors and contains exactly n roots, which 
are not necessarily distinct. Whether these points are real, imaginary or complex, their 


coexistence may create surprising behaviors. 


Every n"-order polynomial possesses exactly n roots 


Any polynomial of the form 


p(x) =a,x" +a, x") +4,_,x"° +...4,x° +a,x+a, 


can always be expressed as 
p(x) = a, (x— 2, (x — 2, Mx— 2,9) X— 2 M*-Z,) 


ma a,] [(-z,) 


where the points z; are the polynomial roots, and they may be real, imaginary or 
complex. | | 


Figure 2.7. Fundamental Theorem of Algebra (After Smith, 1977) 
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With each additional fixed point, coexisting attractors can exhibit varying 
behaviors. Such behaviors are actually emerging in a sort of competitive state — with 
each vying to influence a solution’s trajectory. As Figure 2.8 suggests, such behavioral 
effects may or may not extend globally. Each attracting region is called a basin of 
attraction — the set of starting points whose trajectories are asymptotic to a bounded 
region. Competition amongst the fixed points, in the effect upon x,, exists near and on 


basin boundaries. Moving the fixed points can create new basins and destroy old ones. 





OS pedo’ 





Figure 2.8. Competing Effects of Multiple Fixed Points 


Basin boundaries can take on infinitely many shapes. And basin boundaries can 
be far more complicated than a simple curve, and in most instances are. Within their 
intricate patterns, commonly referred to as Julia sets, is the key that reveals some erratic 
behaviors. Where the basins interact and compete, the behavior is not so obvious. Figure 
2.9 demonstrates how nearby starting points, which are expected to have similar 
behaviors, can assume distinct solution paths, particularly near basin boundaries. 
Consider each starting point, x;, x2, and x3. Despite their ‘nearness’, starting points x, and 


x3 reach a root (through different roots) while x2, which begins on a basin boundary, 


10 





never settles. Hence, behaviors have a sensitive dependence on starting points, and can, 


at times, be considered chaotic. 
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Figure 2.9. Behavioral Effects of Multiple Fixed Points (After Pergler, 1999) 


D. JULIA SETS 

Further consideration of such behaviors begins with observing the role of Julia 
sets. Julia sets are the boundary of basins of attraction -- distinguishing which starting 
points are ‘prisoners’ to some fixed points’ basin, and others that ‘escape’ them. 
Consider the following example in Figure 2.10 (After Peitgen, 1992). Note that 
‘prisoner’ points converge to some basin, while ‘escapees’, had any existed, would never 
settle. While the Julia set may be quite complicated, its role remains crucial in revealing 


the coexistence and competition of complex behavior. 
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Rules: Within the bounded region, select an arbitrary starting point. Move about to 
the next point as indicated (by Newton-Rhapson’s method for cube roots of unity), and 
continue until one, the path halts — in which the next destination is itself (marked by 
X), or two, the path becomes cyclical. Evaluate each starting point. 
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3 4 5 6 7 8 9 1 


> DWOawoeomnmairtrrear er 


faa 
No 


While it is apparent that three basins of attraction exist, there is another valuable 
piece of information. Through adding the number of moves necessary from an 
arbitrary starting point to a fixed point and coloring the basins of attraction, the Julia 
set (approximated by the white boundary) becomes apparent. 





Figure 2.10. Julia Set Example 
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Il. NUMERICAL METHODS 


Since numerical methods are capable of approximating the zeros of an analytic 
function, root-finding methods serve as the tool of choice. Such methods come in many 
shapes and sizes. Some are rather simple; others are complex. Each, however, employs a 
different iterative approach that affects the geometry of numerical solutions, ad thus 
impacts on dynamical interpretations. Consequently, an investigation of such methods 1s 
necessary. 

Although there are many root-finding methods to choose from, their conceptual 
origin is the same — that is, all stem from a successive point-wise approximation of an 
arbitrary function’s root. For example, Taylor’s theorem approximates a function value 
and, when truncated accordingly, is a numerical method of some sort. Figure 3.1 
illustrates another, and more recent, method yielding a one parameter family, e, of single- 
point numerical methods capable of finding roots (After Popovski, 1979). Of particular 
interest were the Newton-Raphson, Chebyshev, Halley, and Laguerre methods. 

For illustrative purposes, a uniform approach is applied. No matter the root- 
finding method, each considers the function, f(z)= z°-1, and is restricted to real arithmetic 
for its geometrical interpretation. From an arbitrary point, x,, approximations are 


computed so as to satisfy Popovski’s single-point method, y(z) = p, + p,(x— p,)°, save 


our final method. Successive computations yield more approximations, xg, x), X2,.... Xn- 
To ensure the dynamic behavior is clear, approximations are represented numerically and 


geometrically. 
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Solving f(z)=0 can be found on the basis of a single-point approximation by the four 
parameter function 


y(Z) = P, + p,(Z- p3)° 
Consider a function f(z) on an interval [a,b] where f(a)f(b)<0 and f'2)f’@A, 
zéfa,b]. Let z;€[a,b] be the i approximation to the root réf[a,b] of f(z)=0. Then the 


following approximation to the root r, 2;+; may be obtained from the system of 
equations where 


¥(Z44;) = 0 
yz, ) = f(z, )= f°, d=, 1, 2. 


and when solved yields 


2444 = 2, + (e-1) 


Pea e ft, se) 7 
fray el fy 


] 

| 

Figure 3.1. Popovski’s Single-Point Iteration Formula 
No matter the function to be approximated, special parameter values, e, reveal 
familiar methods. When e approaches one, the single point iteration formula reduces to 
the popular Newton-Raphson method (Figure 3.2). The approximation method simply 
computes a tangent line to the point x, of our function. When e equals one-half, the 
single point iteration formula reduces to Chebyshev’s method (Figure 3.3). The 
approximation method, rather than line, computes a tangent horizontal parabola to the 
point x, of our function. When e equals negative one, the single point iteration formula 
reduces to Halley’s method (Figure 3.4). The approximation method computes a tangent 
hyperbola to the point x, of our function. Laguerre’s method takes a different approach 
(Figure 3.5)(Press, 1988). Rather than computing a tangent near the point x,, the method 


mimics the function’s behavior there — that is, an 7 order function receives an n" order 
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polynomial approximation. For each of these approximations, the root of its tangent or 


n'® approximation typically represents a better approximation to our function’s root. 


Newton-Rhapson Single Point Approximation 


y(x) = p, + P(x — P3) 
where p,=f(%)— FH %) 
yy f(%, ) 
p; =0 


yielding y(x)= f(x, )+ f(x, (x-x,) 


y(x) approximations... 


n 
0 
1 
2 
ae 
4 
5 
6 





Newton-Rhapson Iterator is 








let F(X) 
*n+] ~*n f"( xy) 
for fix) =x -1... 
G7-] 
Xng] = Xp ~ . 
3k 2 





Figure 3.2. Newton-Rhapson Method 
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Chebyshev Single Point Approximation 


y 
L¢’(x,)P 


where p, = f(x,) +> 


y(x) = p, + p,(x- py 





f(x) 
Pr = 2f (&, a, ~ Ps” 
3 =X, + L(%,) 
2f (x,) 
yielding y(x)= f(x j Ln ear j- LG) : x-x.- f(x, ) 2 
eee ea) 


y(x) approximations... 





_ Chebyshev Iterator is 
_, Lin) fn P fn) 


x = , 
ntl n f (i) 2f'( xp y 
jor fa) =x 1 
3... 3_7)2 
Xia} — , Agel (xp-ly 
3x 2 Ox b 





Figure 3.3. Chebyshev Method 


16 





Halley Single Point Approximation 








P2 
= Pp, + 
y(x) = p Es, 
, P4 
where P, = F(™,) oe 
__4lfr’@v 
= GOT 
came x, + atl On) 
fF (%,) 
, 2 , 3 
yielding y(x)= f(x, ar nae Pa, 
n f(x, pf xn, 2d 
f(x, ) 


y(x) approximations... 





Halley Iterator is 


ex PI OWS Cn) 
mm" OF (xP — FFP O,) 


for fx) =x - 1... 


: 6(x,,"- i ie 
18 x,7 - 6(x,°-1)x,, 


n+] —~ “yn 





Figure 3.4. Halley Method 
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Examining relations between the polynomial, its roots, and its derivatives 




















— io P (x) =(x-x,)(x—-x,)...(x-x,) 
potynomiat. In P,(x)| = In|(x- x) +In|(x- X,) +... + In|(x— x,)| 
d|P(x)| 1 1 1  P 
. ie id ee ee On =—_=G 
Obtain ax X—-X, X—-X, X—-X, P 
derivative ; 2 
relations: d P,(x)| = | | | E P 
~—te = > t¢—— rt... $+; =] | -— = 
dx (x—x,)? (x-x,) (x-x,) [P| P, 


Assume root x, is distance a from current guess, while all other roots are 
assumed to be located at a distance b, where 








a =XxX—-X, b=x-x,, f= 2 5.04n 
1 n-l 
a b Oe auarr: 
a b 


eee =e 
G+,{(n—-1)(nH -G?) p- Fe 


yielding y(x) = [x a (x, —a, ) [x a (x, _ b, yl" 


y(x) approximations... 





Laguerre Iterator is 


Nise ge ore ee 
fl), \( LCD) _(£ODY | (£0) -S 
ro aha LU) Gaal hea " | 


Figure 3.5. Laguerre Method 
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While each method can determine the appropriate root, certain methods are 
preferred. As Figure 3.6 illustrates, when monotonic behavior exists, the preferential 
order is clear — Laguerre, Halley, Chebyshev, and then Newton-Rhapson. In Laguerre, 
Halley and Chebyshev, the approximating curves echo the shape of our function. 
Newton-Rhapson’s approximating curve is restricted to a simple line. Although 
convergence is guaranteed, it varies according to the step sizes of the approximating 
methods. With the smallest step size, Newton-Rhapson is only quadratically convergent 
while the other methods having larger step sizes are cubically scaveiseitt When a 
function is not monotonic, the preference is generally uncertain — with no single method 
consistently better than the others. Clues to determining such an ordering begins with the 


careful observation of each method’s geometry. 


x-] 
Newton 
Chebyshev 
Halley 
Laguerre 





Figure 3.6. Monotonic Behavior & The Methods 
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IV. NUMERICAL METHODS’ GEOMETRY 


Again, numerical methods can sort enous a dynamical system’s behavior 
through finding its roots and examining their affect. With different methods yielding 
different behaviors, an examination of individual numerical methods 1s necessary. Recall 
how our four numerical methods, while obtaining the appropriate root, all sought distinct 
solution paths to it. Consider now a mesh of complex starting points, rather than a single 
real point, with an assortment of fixed points. What are the numerical solution paths 
now? What is the effect of the competition and coexistence of fixed points? What 
numerical method is preferred, if any? Answers to these questions appear when mapping 
and coloring each numerical method’s solutions. 

While an n” order complex polynomial with distinct roots partitions the complex 
plane into n number of basins, the partitions may or may not be equally distributed — or 
even connected for that matter. In an ideal setting, these attracting regions resemble a 
Voronoi diagram — regions containing all points that are the nearest neighbors to the 
polynomial’s zero (Figure 4.1). Few things, though, are ideal. Rather, an attracting 


region contains all starting points that asymptotically approach the zero, despite their 


locality. 
Single oe Multiple F ixed Points 
Fixed Point Ideal Basins (Voronoi Diagrams) Calculated Basins 
e @ ° 4 
° o ek e 


Gia alk Vie aot Seeley wee ee Eee és 
ae 
F 
® @ e 


Figure 4.1. Basins of Attraction 
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One popular method to visualize these regions is basin coloring. The process 
simply assigns n colors to the n basins, executes some numerical method to calculate 
which initial points within a bounded region or mesh converge to a particular basin, and 
paints that basin’s color to that point (Figure 4.2). Furthermore, the number of iterations 
necessary to converge to a root can be shown through variety of color intensities. Points 
calling for fewer iterations appear with greater intensity. Through employing these tools, 


sensible geometric interpretations are possible for nearly all complex polynomials. 


Single Multiple Fixed Points 


Fixed Point Ideal Basins (Voronoi Diagrams) | Calculated Basins 





Figure 4.2. Basins of Attraction Coloring 


A. PURE REAL AND PURE IMAGINARY ROOTS 

When considering pure real or pure imaginary roots, the geometries, while still 
creating a variety of basin shapes, sizes, and complexities, are remarkably similar — only 
differing by a rotation to the appropriate coordinate the axes. Consequently, observations 
for one case support the other. 

Even in what appears to be the simplest of settings, real roots, our basins of 
attraction are not ideal (Figure 4.3). Whether considering the case of equally or unequally 
distributed roots (Figures 4.4 — 4.9), nearest neighbor convergence fails to hold. With the 


exception of Laguerre’s method, the shape of the basins roughly appears to conform to 


Zz 











that of a hyperbola. Laguerre’s basin, on the other hand, resembles our anticipated basin 
shape to a minor degree — especially for those eager to see some relationship. Symmetry 
is another common factor that plays a role in shaping the basins. Equally distributed 


roots generate symmetry throughout the geometry; unequally distributed roots do not. 














Figure 4.3. Ideal Basins & Associated Roots 


With slight exceptions sone basin boundaries, the basin sizes for these are fairly 
comparable. For each, there exists some effective radius of convergence — that is, points 
in the neighborhood of a root tend toward that particular root (Figure 4.4). As Figures 
4.4 — 4.9 suggest, that effective convergence radius not only changes amongst the 
numerical methods applied, but also with each polynomial considered. With higher order 
polynomials commonly creating more and more complex geometries, such radii are often 
greatly reduced. 

Each method also bears some sensitive dependence on starting conditions. With 
nearby starting points assuming distinct solution paths, unpredictable behaviors can 
result. Although Figure 2.7 pointed out the concept initially, chaos’ impact cannot be 
ignored — particularly with it present in every method’s geometry. Figure 4.10 further 
reveals that basin boundaries may be self-similar, with infinite levels of detail, 


characteristic of fractal geometry. 
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Figure 4.4. Equally Spaced Roots — 3" Order 
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Figure 4.5. Equally Spaced Roots — 4" Order 
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Figure 4.6. Equally Spaced Roots — 5" Order 
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Figure 4.7. Unequally Spaced Roots — 3 Order 
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Figure 4.8. Unequally Spaced Roots ~ 4" Order 
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Figure 4.9. Unequally Spaced Roots — 5" Order 
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Figure 4.10. Chaos Everywhere 


In considering the sensitive dependence on starting conditions, one need to only 
observe the ‘decorations’ along the basin boundaries for each method’s geometry in 
terms frequency, size, and structure. As a consequence of the competition and 
coexistence of more and more fixed points, the decorations appear with greater frequency 
with higher order polynomials, yet their size decreases. Whether equally or unequally 
spaced real roots, the Newton-Rhapson, Chebyshev and Halley methods appear to 


experience chaotic dynamics in a phase-shifted manner with each other — an unexpected 
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outcome of their iterators. With rather clean, crisp boundaries, Laguerre’s approximation 

technique provides better, though not absolute, predictability over the other methods. 
While pure roots create nifty geometries, things get really interesting with mixed 

roots. 

B. MIXED ROOTS 


Mixed roots, those roots containing both real and imaginary components, provide 
a rich variety of basin shapes, sizes, and complexities. In many instances, there are 
striking similarities amongst these types of roots and pure ones. But when differences 
appear, a spectrum of spectacular geometries develops. 

1. Roots of Unity (Equally Distributed Roots) 

In the simplest of settings, roots of unity, there are more similarities than 
differences. Nearest neighbor convergence fails to hold, save Laguerre’s approximation 
to the third order polynomial (Figure 4.11). As expected from the equal distribution of 
roots, basin shapes are symmetric (Figures 4.11 — 4.13). Again, these shapes lend to the 
equally distributed sizes of each basin. 

Basin boundaries vary considerably — spanning from the very simple to the 
exceptionally intricate. Consequently, basins are more disconnected. As for sensitive 
dependence on starting conditions, the ‘decorations’ for each method’s geometry in terms 
frequency, size, and structure remains similar to the case of real roots. Furthermore, the 
effective radius of convergence is also affected accordingly (Figure 4.11). Figures 4.11 - 
4.13 suggest, that effective convergence radius not only changes amongst the numerical 


methods applied, but also with each polynomial considered. 


31 











Again, as larger order polynomials are considered, basins become more 
complicated. Such behavior occurred previously, so this is of no great surprise. What is 
of great surprise is the rapid degradation of the Laguerre method near the origin (Figure 
4.14). The apparent ‘disks of chaos’ seem analogous to Feigenbaum’s universality — that 
qualitative changes leading from order to chaos and chaos into order exist (Peitgen, 


1992). 
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Figure 4.11. Roots of Unity — 3" Order 
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Figure 4.12. Roots of Unity — 5" Order 
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Figure 4.13. Roots of Unity — 7" Order 
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Figure 4.14. Roots of Unity — 7" Order (Zoom) 
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2s Unequally Distributed Roots 

When roots are positioned irregularly, interesting things can and do happen. In 
general, the poneeete of nearest neighbor convergence failing, basin ape size, and 
complexity are as noted previously — but perhaps in a more pronounced fashion. 

Third order polynomials produce a variety of familiar behaviors. While Figures 
4.15 & 4.16 echo the behaviors previously found with real roots, Figures 4.17 & 4.18 
behave similarly to the roots of unity. Why the difference? Simply put, one is nothing 
more than a skewed version of the other. When one fixed point is a ‘near-enough’ 
reflection of another, a ‘near-enough’ symmetric geometry results; otherwise, a distorted 
geometry develops. Preference for a particular numerical method is subject to the 
presence of such symmetry. 

Figures 4.19 — 4.24 reveal how convergence changes with the next order of 
polynomials. With Figure 4.19 roots lying on a straight line, it is nothing more than a 
rotation of Figure 4.5, and it assumes similar behaviors. Figures 4.20 — 4.22, while 
containing a reflection of a fixed point to another, contain irregular and unexpected 
basins shapes and decorations — resulting from the competition and coexistence of a 
fourth fixed point. Halley’s method generates the most pronounced irregularities, to 
include distinct breaks in basin connectivity. Figures 4.23 & 4.24 yield expected 
behaviors, save Halley. 

As for hi gher order polynomials, Figures 4.25 — 4.33 depict various geometries 
for fifth and sixth order polynomials. While the geometries are qualitatively different, 
their interpretation is found through the application of previous geometry’s behaviors 
(Figures 4.4 — 4.24). 
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In all these instances, basin shapes, sizes and complexities vary considerably. 
And in nearly every case, the geometry remains unpredictable. But within this chaotic 


environment, is there any order? 
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Figure 4.15. Mixed Roots — 3 Order 
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Figure 4.16. Mixed Roots — 3 Order 
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Figure 4.17. Mixed Roots — 3™ Order 
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Figure 4.18. Mixed Roots — 3™ Order 
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Figure 4.19. Mixed Roots — 4" Order 
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Figure 4.20. Mixed Roots — 4" Order 
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Figure 4.21. Mixed Roots — 4" Order 
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Figure 4.22. Mixed Roots — 4" Order 
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Figure 4.23. Mixed Roots — 4"" Order 


45 











Newton-Rhapson 





Chebyshev Laguerre 





Roots 


9, i, -6-6i, —5+4i 








Figure 4.24. Mixed Roots — 4" Order 
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Figure 4.25. 5" Order Mixed Roots 
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Figure 4.26. Mixed Roots — 5" Order 
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Figure 4.27. Mixed Roots — 5" 
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Figure 4.28. Mixed Roots — 5" Order 
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Figure 4.29. Mixed Roots — 5" Order 
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Figure 4.30. Mixed Roots — 6" Order 
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Figure 4.31. Mixed Roots — 6" Order 
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Figure 4.32. Mixed Roots — 6" Order 
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Figure 4.33. Mixed Roots — 6" Order 
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V. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

The Newton-Raphson, Chebyshev, Halley and Laguerre approximation methods, 
serve as powerful tools in evaluating complex polynomials’ roots. These different 
methods, however, can yield different solutions from identical starting points. In 
determining any preference for the numerical methods, consideration must be given to 
the polynomial at hand, when do root finding methods converge and how long for 
convergence. 

Whether low or high order, the Laguerre approximation method tends to fare 
better than other methods. With relatively simple basin boundaries, the method not only 
affords a greater effective radius of convergence but also increased behavior 
predictability. In many instances, the Newton Raphvon and Halley geometries are nearly 
indistinguishable for the same reason. When fixed points are a reflection of another, 
Halley’s method assumes a much larger effective radius over the Newton-Raphson 
method, and it can even outdo Laguerre’s method. Chebyshev’s method, filled with 
complex boundaries and relatively small effective radii, remains the worst of the group. 

With the methods relatively comparable 1n computational speed, the greater 
emphasis rests in basin shape, size, and complexity. 

Ultimately though, the method of choice depends on the complex polynomial at 


hand. 
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B. RECOMMENDATIONS 
While room for further research in this topic exists, a particular effort with respect 
‘to more numerical methods, calculating effective radii, and consideration for repeated 


roots would prove both challenging and rewarding. 
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APPENDIX. BASIN CODE (MATLAB) 


6 FE 8 A He HE EE Ee ee ee ee HC He ee C2 Ce HE HE 2 2 eK OR EE 2 A EE 2 IC 2 C2 C2 He a 2 7, 
% This MATLAB program computes basins of attractions for complex, analytic 
% polynomials using various numerical methods. These methods include 

% Newton-Rhapson, Chebyshev, Halley, and Laguerre. 


% 
Jo 
% 
No 
% 
%o 
% 
Yo 
% 
% 
% 
% 
Jo 
% 
% 
% 
% 
%o 
Jo 
% 
Jo 
Jo 
% 
% 
%o 


User inputs: f: Analytic function: 


Notes: 


ie. f=[1001] ==> 243 +1 

method1: Numerical Method 

i.e. N’= Newton-Rhapson, 'C’ = Chebyshev 

H’ = Halley, L’ = Laguerre 
t_limit: Maximum acceptable absolute difference between the 
computed and actual root for both axis. 

i.e. tol=.01+.01i 

max_iteration: Maximum number of iterations before starting 
point becomes a member of the Julia set 
i.e. max_iteration=100 


With an nth degree polynomial generating n roots, the user must 
include n sets of the following codes to account for all roots. 


if abs(p_n-actual_root(2)) <= tol 
break 
end; 


if abs(p_n-actual_root(1)) <= tol 
root_color_code(real_counter,imag_counter)=1; 
end; 


Also, nth degree polynomial requires n+1 color assignments 


% 
% 
% 
Jo 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 


To 76 76 A 2 HR 2 AE OE EE EE He 2 ae 2 AG 2 ae A ee eR ee Ee Ee Ie Ee IE Ee ee I ee 2 eK AK OP, 


function basin_generator=basins(f,method1,t_limit,max_iteration) 


% Defining/reseting initial conditions 


iteration=0; 


iteration_counter=0; 


roots_found=0; 


imag_counter=0; 


real_counter=0; 


n=length(f)-1; 


method=char(method1); 
tol=t_limit +t_limit*i; 
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% Determining the actual roots, coloring assignments, and complex plane 

% boundaries and starting point step size 

actual_root=roots(f); 

root_colors=([{1,0,0; 0,1,0; 0,0,1; 1,1,0; 1,0,1; 0,1,1; .5,1,0; 1,.5,0; .5,.5,.5; 1,1,1]); 

bound=max([abs(max(imag(roots(f)))) abs(min(imag(roots(f)))) abs(max(real(roots(f)))) 
abs(min(real(roots(f))))]); 

imag start_pt=-bound-.5; 

imag end _pt=bound+.5; 

real start _pt=-bound-.5; 

real end_pt=bound+.5; 

step=bound/30 


% Imaginary axis boundaries/do-loop 

for imag _axis=imag_start_pt:step:imag end pt 
imag counter=imag counter+] 
real_counter=0; 


% Real axis boundaries/do-loop/assigning starting points 
for real_axis=real_start_pt:step:real_end_pt 

real _counter=real_counter+1; 

p_n_ 1=real_axistimag axis*1; 


% Resetting iteration counter/root 'a' measure for next starting point 
iteration=0; 
a=tol+1; 


% Iteration to determine convergence or Julia set member 
while iteration<=max_ iteration 


% Fail safe -- No iteration necessary is starting on a root 


if polyval(polyder(f),p n_1)==0 
break 
end 


% Newton-Rhapson Iterator 


if method='N' 
p_n=p_n_1-polyval(fp_n_1)/polyval(polyder(f),p n_1); 
end 


% Chebyshev Iterator 
if method='C' 
p_n=p_n_1-polyval(f,p_n_1)/polyval(polyder(f),p_n_1)- 
((((polyval(fp_n_1))*2)*(polyval(polyder(polyder(f)),p_n_1))) 
/(2*(polyval(polyder(f),p_n_1))*3)); 
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end 


% Halley Iterator 
if method="H' | 
numerator1=2*polyval(f£p_ n_1)*polyval(polyder(f),p n_1); 
denominator1=2*(polyval(polyder(f),p_n_1))*(polyval(polyder(f),p_ n_1))- 
(polyval(f;p _n_1))*polyval(polyder(polyder(f)),p_n_.1); 
p_n=p_n_1-numerator1/denominator1; | 
end 


% Laguerre Iterator 
if method="'L' 
if iteration ~= 0 
_n_1=p_n; 
else 
p_n_1=real_axist+(imag_axis)*1; 
end 


% Fail safe -- No iteration necessary is starting on a root 
if polyval(f;p n_1)==0 
break 
elseif polyval(f,p_n_1)~=0 
G=(polyval(polyder(f),p_n_1)/polyval(f,p n_1)); | 
H=G“2- polyval(polyder(polyder(f)),p_n_1)/(polyval(fp_n_1)); 
if (G + sqrt((n-1)*(n*H-G“2)))==0 
break 
end 
if (G - sqrt((n-1)*(n*H-G“%2)))==0 
break 
end 
if abs(G + sqrt((n-1)*(n*H-G%2))) > abs(G - sqrt((n-1)*(n*H-G“2))) 
a=n/(G + sqrt((n-1)*(n*H-G‘2))); 
else 
a=n/(G - sqrt((n-1)*(n*H-G‘2))); 
end 
end 
end 


% Updating computed roots/iteration 
iteration=iteration+1; 
if method=="L' 
p n=p_n 1-a; 
else 
pn _1=p n; 
end 
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% Is computed root within tolerance of an actual root? 
if abs(p_n-actual_root(1)) <= tol 
break | 
end 
if abs(p_n-actual_root(2)) <= tol 
break 
end 
if abs(p_n-actual_root(3)) <= tol 
break 
end 


% Extra statements for n roots 

% if abs(p_n-actual_root(4)) <= tol 
% break 

% end 

% if abs(p_n-actual_root(5)) <= tol 
% break 

% end 

% if abs(p_n-actual_root(6)) < tol 
% break 
% end 

% if abs(p_n-actual_root(7)) <= tol 
% break 

% end 

% if abs(p_n-actual_root(8)) <= tol 
% break 

% end 

% if abs(p_n-actual root(9)) <= tol 
% break 

% end 


% No root found, update variable counters for next iteration 
end % while 


% Iteration/Computed root trackers 
iteration _counter(real_counter,imag_counter)=iteration; 
computed_root(real_counter,imag_counter)=p_n; 


% If computer root within tolerance of an actual root, do color assignment? 
if abs(p_n-actual_root(1)) <= tol 
root_color_code(real_counter,imag_counter)=1; 
elseif abs(p_n-actual_root(2)) <= tol 
root _color_code(real_counter,imag_counter)=2; 
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elseif abs(p_n-actual_root(3)) < tol 
root_color_code(real_counter,imag_counter)=3; 

% Extra statements for n roots 

% elseif abs(p_n-actual_root(4)) <= tol 

%  root_color_code(real_counter,imag counter)=4; 

% elseif abs(p_n-actual root(5)) < tol 

% root color code(real_counter,imag_counter)=5; 

% elseif abs(p_n-actual_root(6)) <= tol 

%  root_color_code(real_counter,imag_counter)=6; 

% elseif abs(p_n-actual_root(7)) <= tol 

%  root_color_code(real_counter,imag_counter)=7; 

% elseif abs(p_n-actual_root(8)) <= tol 

% root_color_code(real_counter,imag_counter)=8; 

% elseif abs(p_n-actual_root(9)) < tol 

%  root_color_code(real_counter,imag_counter)=9; 

else 
root_color_code(real_counter,imag_counter)=10; 

end % if 


end % forreal axis 
end % forimag axis 


% Building the true color specification for root_color_code using 
% anm-by-n-by-3 array of RGB values. 

b(:,:,1) = iteration_counter; 

b(:,:,2) = 1teration_counter; 

b(:,:,3) = iteration_counter; 


% Scaling the colors to include iteration levels 
for 1= 1:length(root_color_code) 
for j = 1:length(root_color_code) 
bGj,1,:)}= root_colors(root_color_code(i,j),:)*(((iteration_counter(i,j) 
/(max(max(iteration_counter)))))*.37); 
end 
end 


% Draw figure with appropriate title 
figure(1) 
image(b) 
image(b,'XData',[-bound-.5 bound+.5],"YData',[-bound-.5 bound+.5]); 
if method='N' 
title("Newton-Rhapson’); 


end 
if method="C' 
title(‘Chebyshev’); 
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end 

if method—"H' 
title("Halley'); 

end — 

if method=="L' 
title(‘Laguerre’); 

end 

axis Square 
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